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Abstract 

We study the steady state properties of a 2D granular mixture in the presence of energy driving 
by employing simple analytical estimates and Direct Simulation Monte Carlo. We adopt two dif- 
ferent driving mechanisms: a) a homogeneous heat bath with friction and b) a vibrating boundary 
(thermal or harmonic) in the presence of gravity. The main findings are: the appearance of two 
different granular temperatures, one for each species; the existence of overpopulated tails in the 
velocity distribution functions and of non trivial spatial correlations indicating the spontaneous 
formation of cluster aggregates. In the case of a fluid subject to gravity and to a vibrating bound- 
ary, both densities and temperatures display non uniform profiles along the direction normal to the 
wall, in particular the temperature profiles are different for the two species while the temperature 
ratio is almost constant with the height. Finally, we obtained the velocity distributions at different 
heights and verified the non gaussianity of the resulting distributions. 

PACS numbers: 02.50.Ey, 05.20.Dd, 81.05.Rm 
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I. INTRODUCTION 



Granular materials present a rich and intriguing phenomenology, which has attracted the 
interest of the scientific community since the nineteenth century [[L]]. However, in spite of its 
recent progress the theoretical study of granular gases, i.e. of fluidized granular particles, is 
certainly less advanced than that concerning ordinary molecular fluids. The reason for this 
state of affairs is the presence of dissipation due to inelastic collisions and of friction with the 
surroundings, which prevents these system to reach thermodynamic equilibrium. In fact, in 
order to render stationary a granular system one needs to inject continuously energy into 
the system. This can be done, for instance by shaking or vibrating the grains. 

In the present paper we illustrate the results of a numerical investigation concerning the 
properties of a two component granular mixture, modeled, following an established tradition, 
by inelastic hard spheres (IHS) with different masses, restitution coefficients, radii, and 
subject to different forms of external drive. The physical motivation for our study stems 
from the fact that in nature most granular materials are polydisperse from the point of view 
of their sizes and/or of their physical and mechanical properties. The theoretical study of 
granular mixtures has attracted so far the attention of several researchers 0, [|, |], [|. These 
studies comprehend both freely cooling and uniformly heated granular mixtures and have 
been performed almost contemporaneously with laboratory experiments |J [7J. The most 
striking outcome is the lack of energy equipartition, i.e. the presence of two different kinetic 
temperatures, one for each species. 

In the present paper our interest will be concentrated on the stationary state obtained 
by applying an energy feeding mechanism, represented by a uniform stochastic driving force 
or by a vibrating wall. In the first case the work of Barrat and Trizac {§] and Garzo et 
al. H provides a simple theoretical account, based on the Boltzmann equation approach 
and on the Direct Simulation Monte Carlo method (DSMC) || |§, of the lack of energy 
equipartition. 

The present work differs from the recently appeared papers dealing with granular mixtures 
in several aspects: 

• the heat bath, which drives the system, is at finite temperature, since we have intro- 
duced a finite friction of the particles with the surroundings and we have considered 
the possibility of having spatially inhomogeneous states; 
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• in addition, we considered a second situation in which the particles subject to a vertical 
gravitational field receive energy anisotropically from the bottom vibrating wall of the 
container; 

• in both cases we provide information about the presence of inhomogeneities in the 
system, e.g. density clustering and non-uniform density and temperature profiles. 

In Sec. II we present the model fluid and two different mechanisms of energy supply. In 
sect III we discuss the first sub-model, the one with the heat bath, and obtain by means of 
an approximate analytic method an estimate of the partial temperature of each component. 
Subsequently we study the same sub-model with a DSMC algorithm. In sect. IV we study 
numerically the second sub-model (the one with gravity and vibrating wall) by means of 
DSMC simulations. Finally in the last section, V, we present our conclusions. 

II. DEFINITION OF THE MODELS 

We shall consider a dilute inelastic gas constituted of Ni particles of mass mi and iV 2 
particles of mass m 2 subject to some kind of external driving (this will be specified in the 
following). We suppose that the interactions between the grains can be described by the 
smooth inelastic hard sphere model (IHS) jnj, thus we specify only the radius of the spheres, 
their masses and the fraction of the kinetic energy dissipated at each collision. This can 
be done by defining three different restitution coefficients a^-, i.e. au,« 2 2, and «i 2 = a 2 i, 
which account for normal dissipation in collisions among particles of type % and j. No 
internal degrees of freedom (e.g. rotations) are included. 

One can describe the velocity changes induced by the instantaneous inelastic collisions of 
smooth disks labeled 1 and 2 of diameter o\ and a 2 by the following equations: 



where n = 2(xi — x 2 )/(a Kl + cr K2 ) is the unit vector along the line of centers xi and 
x 2 of the colliding disks at contact and Ki,k 2 are the species (1 or 2) to whom particles 




(lb) 



(la) 
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1 and 2 belong. An elementary collision conserves the total momentum and reduces the 
relative kinetic energy by an amount proportional to (1 — ck^ s )/4. The collision rule we 
have adopted excludes the presence of tangential forces, and hence the rotational degrees of 
freedom do not contribute to the description of the dynamics. 

Since the particles suffer mutual collisions and loose kinetic energy, in order to achieve 
a steady state, one needs to supply from the exterior some energy. The energy source has 
been modeled in two different fashions. 

In the first sub-model we have assumed that the particles, besides suffering mutual col- 
lisions experience a uniform stochastic force and a viscous damping [|TT1 , |T2|1 . The presence 
of the frictional, velocity-dependent term in addition to the random forcing, not only is 
motivated by the idea of preventing the energy of a driven elastic system (a K1K2 — > 1), to 
increase indefinitely, but also mimics the presence of friction of the particles with the con- 
tainer. Moreover a fluctuation dissipation relation is assumed between the viscous force and 
the intensity of the noise. Even in extended systems with small inelasticity the absence of 
friction may cause some problems of stability . We shall be mainly interested in the sta- 
tionary situation determined by the balance between the energy feeding mechanism and the 
dissipative forces due both to the inelasticity of the collisions and to the velocity-dependent 
frictional force. 

In the second sub-model the grains are constrained to move on a frictionless inclined 
plane and the bottom boundary vibrates (as a thermal [I3| or deterministic [14" oscillating 
wall) giving energy only to the particles bouncing on it. Periodic boundary conditions are 
assumed laterally. 

Since we consider throughout only sufficiently low density systems successive binary col- 
lisions are effectively uncorrelated and Boltzmann equation can be used to describe the 
non equilibrium dynamics. In fact, by assuming the validity of the Boltzmann molecular 
chaos hypothesis, introduced to treat the collision term, it is straightforward to derive the 
governing equations for the probability density distribution function of each species. More- 
over in this situation we are allowed to perform Direct Simulation Monte Carlo to integrate 
numerically the inhomogeneous Boltzmann equation of the sub-models. 
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III. UNIFORMLY HEATED SYSTEM 



In order to see the effect of the heat bath let us consider the system in the absence of 
collisions. In this case, the evolution of the velocity of each particle is described by an 
Ornstein-Uhlenbeck process. If we require that the two components must reach the same 
granular temperature in the limit of vanishing inelasticity we have two different possibilities 
to fix the heat bath parameters: 



d t x i (t)=v i (t) (2) 



rriidtVi (t) = - 7 v, (t) + y f 2^f b & (t) (3a) 

mid t Vi(t) = -miPvi(t) + ^/2m$T h (3b) 

where % = 1,2 and T), is the heat bath temperature and £(£) is a Gaussian noise with the 
following properties: 

< m > = (4a) 

< 6(*i)&(*2) > = S(ti - t 2 )8ij (4b) 

The associated Fokker-Planck equations for the two cases are respectively: 

dtfi(T, v, t) = ^-V v (yfi(T, v, t)) + ^V^,(r, v, t) + W r / 4 (r, v, t) (5a) 
rrii m\ 

dtfi{T, v, t) = /3V v (vfi(T, v, t)) + ^ V2/i(r, v, t) + vV r /i(r, v, t) (5b) 



A. Spatially Uniform solutions 

When we take into account collisions among particles equations (HP become two coupled 
Boltzmann equations modified by the presence of a diffusion term due to the thermal noise. 
We shall first consider a spatially homogeneous case in order to derive the temperature 
of each species in the homogeneous stationary state. This represents a sort of mean field 
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approximation. In fact, in the original Boltzmann equations the collisions are considered 
to occur only between particles at the same spatial location, whereas here these can occur 
between arbitrary pairs of particles regardless their spatial separation. In the present sub- 
section we compute the second moments, < v\ >, of the distribution functions in order to 
determine the partial temperatures and their ratio. Although the method of derivation of 
the equations for the partial temperatures is not original, we present it in order to render 
the paper self-contained and because it shows the differences between the particular heat 
bath we employed and those chosen by other authors @]. First, indicating by rij = Ni/V 
the partial density of species i, we notice that both eqs. © possess the same equilibrium 
solution: 

although the relaxation properties are different. Only upon adding the inelastic collision 
term the two species display different temperatures. The resulting Boltzmann equation for 
a granular mixture is || ||, || : 

dtfiin; t) = J2 ^-M/i, fj] + tfv'jt + AV„ ■ (vx/0 (7) 

3 

where we have used a compact notation to represent the two different choices of heat bath: 
• Case 1 



f 2 



27^ 
■1 



Pi 



1 

in, 



(8) 



• Case 2 



SOi 



2/3T b 



(9) 
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and Jij[vi\fi, fj] is the collision integral: 



JyM/i, fi\ = 4 I dv 2 I daQ(a • g 12 )(a • s a )[atff i {v , 1 )f j {y' a ) - /,i v, v 2 ) (10) 

The primed velocities are pre-collisional states, which can be obtained by inverting eqs. ([[]) 
Due to the presence of the heat bath terms the system reaches asymptotically a steady 
state, characterized by time independent pdf's. By imposing the vanishing of the time 
derivatives we obtain: 



E JyW/o U\ + Y v " fi + AV " ' (vi/i) = 

3 

or after integrating over v x : 

e y ^ 1Wl 2 ^[vxi/i, /,] + ^ y dvi^v^/i + a y . ( V1 /o = o (12) 

After simplifying the second and the third integral by integration by parts and using the 
normalization property j fidvi = rij we find: 



E / dv 1 vlJ ij [v 1 \f i , fj] +n<^ - 2pi / dv l v 2 l f i {y 1 
j J J 

By defining the partial temperature: 



riiTi = - j d\ 1 m i vlf i 



we can rewrite eq. (|T^) in terms of the temperatures : 



rrti 



2dPi \ n. 



~Y1 / dv l V l J iA V l\fii fj\ + d & 
"H J 



(13) 



(14) 



(15) 



Eq. (|i~5"|) determines the partial temperatures once the fi are known. In practice one can 
obtain an estimate of Tj by substituting two Maxwell distributions: 



d 

(Tfl ' \ ^ 
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and performing the remaining integrals (see @, ||) one gets: 



dT(d/2) 



2/3 i (T 6 -T i ) =<r i n 



2(1 - a? 



m 

/2Ti 27! 



3/2 



3/2 



\ rrii rrij 



mi + rri2 



2% 2T„ 



+ 



nii 



rriA 



1/2 



(16) 



where \iij = rrii/ (m, + rrij). 

By solving numerically the nonlinear system of eqs. flTE| ) one obtains the steady values 
of the partial temperatures in the spatially homogeneous situation. 

In figs. [I] and ^| we report the temperature ratio ^ as a function of a common restitution 



coefficient a, having chosen equal coefficients an 



"22 



ai2 = a. Assuming identical 



concentrations, and a varying the mass ratio we studied the cases 1 and 2. 



— m^O.5, m 2 =l 
■ ■■■ m^l, m 2 =3 
-- m=3, m 2 =l 




FIG. 1: Granular temperature ratio T1/T2 vs. a obtained with the heat bath, case 1, homogeneous, 
using Tb = 1, 7 = 0.1 and different mass ratios. 

We notice that with the second recipe (case 2) the temperature ratio is an increasing 
function of the mass ratio niijm^. this is exactly the opposite of what happens for the case 
1. The experimental observation suggest that the trend of case 2 is the right one. In 
other words, in the case 2 both the friction term and the power supplied are proportional 
to the mass of the two species and such a property is probably true also in the experiments. 
This is why we shall use the recipe of case 2 in the following DSMC simulations. 
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FIG. 2: Granular temperature ratio T1/T2 vs. a using the second recipe for the heat bath, case 
2, homogeneous, with Tf, = 1, (3 = 0.1 and various mass ratios. 

B. Non uniform solutions 

In the present section we illustrate the results obtained simulating the system with the 
heat bath (with recipe 2, i.e. equation (|3l5D) by the so called Direct Simulation Monte Carlo 
according to the implementation described in |TEfl. 



At every time step of length At each particle is selected to collide with a probability 
p c = At/r c (where r c is an a priori fixed mean free time established consistently with 
the mean free path and mean squared velocity) and seeks its collision partner among the 
other particles in a neighborhood of radius tb, choosing it randomly with a probability 
proportional to their relative velocities. Moreover in this approximation the diameter o is 
no more explicitly relevant, but it is directly related to the choices of p c and r# in a non 
trivial way: in fact the Bird algorithm allows the particles to pass through each other, so 
that a rigorous diameter cannot be defined or simply estimated as a function of p c and r#. 
In this section, to indicate the degree of damping, we give the time Tj, = 1/(3 instead of 
(3. This is useful to appreciate the ratio between the mean collision time r c and the mean 
relaxation time due to the bath which indeed is t^. 

In the present section we chose Ni = N 2 = 500 and T b — 1, and equal restitution 
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m 2 /m i 



FIG. 3: Granular temperature ratios T1/T2 vs. mass ratio milml for a binary mixture (DSMC 
simulation) with different values of a, N± = N 2 = 500, L 2 = 1000, If, = 1, Tf, = 10, r c = 0.16, and 
recipe 2 
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FIG. 4: Granular temperature ratios T\jTi vs. restitution coefficient a for a binary mixture 
(DSMC simulation) with different values of m 1 /m 2 , N x = N 2 = 500, I? = 1000, T b = 1, r 6 = 10, 
r c = 0.16, and recipe 2 
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coefficients for all collisions and r c = 0.16. As illustrated in figures |3| and [|the two compo- 
nents display different granular temperatures in agreement with the analytical predictions 
of the homogeneous Boltzmann equations. We also comment that the ratio Ti/T 2 is very 
sensitive to the mass ratio and much less to differences in the restitution coefficients and 
concentrations. This is in agreement with experimental observations |J. 

At a finer level of description we consider the rescaled velocity pdf for different values of 
the inelasticity parameter and different mass ratios 7712/ Vni = 2. One sees that not only the 
deviations from the Gaussian shape become more and more pronounced as we increase the 
inelasticity parameter (1 — a), but also the shape of the two distributions differ appreciably 
in the tails even after velocity rescaling to make the two pdf 's have the same variance. Such 
a property is similar to the one already reported in ||. One can also observe that the 
rescaled pdf of the lighter species has slightly broader tails. The mass ratio also controls the 
deviations from the Gaussianity of the velocity pdf 's. Whereas it is known that the departure 
from the Maxwell-Boltzmann statistics is triggered by the inelasticity of the collisions, the 
larger the inelasticity the stronger the deviation, to the best of our knowledge this is the 
first numerical evidence of the phenomenon, predicted within a Maxwell model in ref. |9. 
Comparing figs. ^ and || one sees that the mass asymmetry enhances the non gaussianity of 
the pdf. 

We have also studied the limits of low and high r&, in figure [7|, to show how the velocity 
distributions change. For values of the characteristic time of the heat bath, 77,, comparable 
with the collision time, r c , the dynamics is essentially controlled by the stochastic accelera- 
tion term. This fact renders the two partial temperatures very close and makes the velocity 
distributions nearly Maxwellian. As r& increases we have observed that the energy dissipa- 
tion due to the inelasticity makes the temperatures of the two species different. Moreover, 
the temperature ratio displays power law decreasing trend as a function of r& whose strength 
depends on the mass ratio (see figure ||). 

In order to obtain information about the spatial structure of the mixture we have per- 
formed an analysis of the following correlation function, already introduced in the context 
of granular media by [fTT , |T2| HJ : 



^(r) = ^3 Ty E ^-|<- x ?l) ( 17 ) 
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FIG. 5: Rescaled (to have variance 1) velocity distributions P(v) vs. v in the numerical experiment 
with the thermal bath (DSMC simulation), for binary mixtures of particles with mass m\ = 0.5 
and m 2 = 1, with N x = N 2 = 500, L 2 = 1000, T b = 1, r b = 10, r c = 0.16, recipe 2, and different 
values of a 

For a spatially homogeneous system we expect that C a p(r) « r d2 with d% = d the dimen- 
sion of the embedding space. However, in the cases we have studied, we found that d 2 < d 
(see fig g), a signature that the system tends to clusterize. The larger the value of t& the 
stronger the deviation. 



IV. SYSTEM WITH GRAVITY AND HEATED FROM BELOW 



We turn, now, to illustrate the results relative to an inelastic mixture subject to gravity 
and confined to a vertical plane of dimensions L x x L y . In the horizontal direction x we 
assumed periodic boundary conditions. Vertically the particles are confined by walls. Energy 
is supplied by the bottom wall vibrating stochastically or periodically according to the 
method employed in |ll| for a one component system. 

The vibration can have either a periodic character (as in ref. [I8|) or a stochastic behavior 
with thermal properties (as in [13[|). In the periodic case, the wall oscillates vertically with 
the law Y w = A w sm(u w t) and the particles collide with it as with a body of infinite mass 
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FIG. 6: Rescaled (to have variance 1) velocity distributions P(v) vs. v in the numerical experiment 
(DSMC simulation) with the thermal bath, for binary mixtures of particles with mass mi = 0.5 
and m 2 = 5, with iVi = N 2 = 500, I? = 1000, T b = 1, r 6 = 10, r c = 0.16, recipe 2, and different 
values of a 

with restitution coefficient «„, so that the vertical component of their velocity after the 
collision is v' = —a w v y + (1 + a w )V w where V w = A w u w cos(u w t) is the velocity of the 
vibrating wall. In the stochastic case we assume that the vibration amplitude is negligible 
and that the particle colliding with the wall have, after the collision, new random velocity 
components v x G (— oo, +oo) and v y G (0, +oo) with the following probability distributions: 



P(v v ) = ^exp(-^r) (18) 

1 W Zl yj 



P(v x ) = ; 1 exp(-^-) (19) 

In this model we assume that the particle do not feel any external (environmental) friction. 
In all simulation we have chosen N x = N 2 — 100, g — 1, mi = m 2 /2 = 1/4, L x = L y = L = 

Vn. 

The Boltzmann equations for the partial distribution functions f K (r,v;t) with k — 1,2 
read 
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FIG. 7: Rescaled (to have variance 1) velocity distributions P(v) vs. v in the numerical experiment 
(DSMC simulation) with the thermal bath, for binary mixtures of particles with mass mi = 0.5 
and m 2 = 5, with JVi = N 2 = 500, L 2 = 1000, a = 0.2, T b = 1, different values for r b , r c = 0.16, 
with recipe 2, and different values of a 



- + v • V r + 9i — ) f K (r, v,t) = Yl J M«> M- ( 2 °) 
«/ ^ 

The simulation results show that the tendency of the grains to form clusters is enhanced 
by such a choice of driving mechanism with respect to the homogeneous heat bath of the 
previous section. In the latter the noise acting uniformly was more effective in breaking the 
clusters. In addition the gravitational force tends to group the particles in the lower portion 
of the container for not too large driving frequencies jnSfl. 



Fig [10] illustrates the partial density profiles and granular temperature profiles in the 
presence of a thermal wall of intensity T p = 20. The density profiles differ slightly near 
the bottom wall where both present a maximum (see inset). The temperature profiles are 
different, putting again in evidence a strong lack of equipartition in the system. Interestingly, 
the temperature ratio is almost constant along the vertical direction, notwithstanding the 
partial profiles are non constant. We have also performed numerical simulations with a 
harmonically vibrating wall: as the maximum velocity of the vibrating wall (which, for 
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FIG. 8: Temperature ratios T1/T2 vs. the viscosity time t& in the numerical experiment (DSMC 
simulation) with the thermal bath for a binary mixture of particles with mass mi = 0.5 and m,2 = 5, 
with Ni = N 2 = 500, L 2 = 1000, a = 0.2, T b = 1, different mass ratios, r c = 0.16, with recipe 2, 
and different values of a 
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FIG. 9: Correlation function used to compute the correlation dimension d,2, C(r) vs. r (see text 
for definition), for a binary mixture (DSMC simulation) with mi = 0.5, m 2 = 5, r c = 0.16, T b = 1, 
Ni = N 2 = 500, I? = 1000, a = 0.2, with different values for r b 
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FIG. 10: Profiles for density and granular temperature n(y), T(y) vs. y for a binary mixture 
with mi = 0.5, rri2 = 2, on an inclined plane with a thermal bottom wall (DSMC simulation) with 
T w = 20, N\ = N2 = 100, L 2 = 200, a = 0.6, a w = 1, g = 1; the dashed lines indicate the ratios 
for both profiles. 

A = 1, is equal to cu) increases, the position of the density maximum raises indicating that 
gravity becomes less and less relevant. As far as the partial temperatures are concerned 
figure (PD illustrates the corresponding situation. 

One sees that the temperatures (see inset) next to the wall attain their largest value, 
then drop to increase again. This indicates that the region far from the bottom is hotter 
because of the lower density of the gas and of the lower collision rate. As the wall temperature 
increases the temperature becomes more homogeneous far from the bottom due to the major 
homogeneity in density. We also notice a small segregation effect. 

Finally, we present the velocity distributions measured in the system with a thermal wall. 
Since the temperature depends on the vertical coordinate y we have computed the velocity 
pdf at different heights. The various pdf are plotted in figure (0) after a suitable rescaling. 
We observe that the distributions deviate appreciably from a Gaussian shape and display 
overpopulated high energy tails. 

Interestingly the value obtained by our simulation for the temperature ratio, having 
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FIG. 11: Profiles for density n(y) vs. y and ratios between granular temperatures of the two 
species T\(y) /T2{y) vs. y, for a binary mixture (DSMC simulation) with mi = 0.5 and m% = 2, 
on an inclined plane with an harmonic oscillating bottom wall with A = 1, different values of u, 
N 1= N 2 = 100, L 2 = 200, a = 0.9, a w = a g = 1. 

chosen — = 0.3 and a single restitution coefficient a = 0.93, 2l ~ 0.75 is not too far from 



the value obtained experimentally by Feitosa and Menon || which is 0.66 ±0.06. One should 
recall that in our "setup" only the lower wall supplies energy to the system, whereas in the 
experiment both walls vibrate. 

V. CONCLUSIONS 

To summarize we have studied the steady state properties of a granular mixture subject 
to two different classes of external drive. 

In the case of an heat bath acting homogeneously on the grains, we have first obtained 
the temperatures of each species by employing an approximate analytic method. We have 
then solved numerically the equations by allowing the density and the temperatures to be 
spatially varying and analyzed the spatial correlations and the velocity distributions. 

Finally, we turned to investigate the properties of the mixture in the presence of gravity 
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FIG. 12: Rescaled distributions (to have variance 1) of velocity P(v) vs. v for particles taken in 
stripes at different heights, for a binary mixture (DSMC simulation) with mi = 0.5, rri2 = 2, on 
an inclined plane with a thermal bottom wall with T w = 20, N% = N2 = 100, I? = 200, a = 0.6, 

= 1. 

and inhomogeneous drive. 

The present results not only confirm the predictions of different temperatures for the two 
species for all types of drive, but also show the existence of different shape functions for the 
velocity distributions and of overpopulated high energy tails in agreement with the findings 
based on Maxwell models ||. 

In order to obtain a direct comparison with experiments we employed parameters com- 
parable to those utilized in the experimental work of ref. || and found temperature ratios 
not too far from the ones there reported. In spite of this success, the agreement seems to 
be limited to the average properties, while it is clear from the observation of the density 
profiles computed numerically that a better treatment of the excluded volume effect is nec- 
essary in order to obtain a more realistic description of the system. Finally, in view of the 
approximations inherent to the Boltzmann approach, we have not tried to include neither 
the effect of rotations of the grains nor the friction with the lateral walls which might be 



both relevant 19 . 
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